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Identifying rare chaotic and regular trajectories in dynamical systems with Lyapunov 

weighted path sampling 

Philipp Geiger and Christoph DellagcQ 

Faculty of Physics, University of Vienna, Boltzmanngasse 5, 1090 Vienna 

Depending on initial conditions, individual finite time trajectories of dynamical systems can have 
very different chaotic properties. Here we present a numerical method to identify trajectories with 
atypical chaoticity, pathways that are either more regular or more chaotic than average. The method 
is based on the definition of an ensemble of trajectories weighted according to their chaoticity, 
the Lyapunov weighted path ensemble. This ensemble of trajectories is sampled using algorithms 
borrowed from transition path sampling, a method originally developed to study rare transitions 
between long-lived states. We demonstrate our approach by applying it to several systems with 
numbers of degrees of freedom ranging from one to several hundred and in all cases the algorithm 
found rare pathways with atypical chaoticity. For a double- well dimer embedded in a solvent, which 
can be viewed as simple model for an isomerizing molecule, rare reactive pathways were found for 
parameters strongly favoring chaotic dynamics. 



I. INTRODUCTION 

The phase space of dynamical systems often exhibit re- 
gions with qualitatively very different dynamics. In the 
Henon-Heiles model or similar low-dimensional Hamilto- 
nian systems, for instance, islands of stability are embed- 
ded in a chaotic sea Other examples for this kind of 
behavior include the Fermi-Pasta-Ulam chain, in which 
special initial conditions lead to physically very differ- 
ent soliton and "chaotic breather" solutions, and gravita- 
tional many-body systems in celestial mechanics, which, 
for appropriate initial conditions, produce orbits that are 
stable for very long times . Although such particularly 
regular (or irregular) trajectories may be very rare, they 
may be responsible for important physical phenomena 
as is the case for chemical reactions, where trajectories 
passing through unstable saddle points regions carry the 
system from one chemical species to another Q • 

Identifying and describing such trajectories is therefore 
of great interest. Recently, Tailleur and Kurchan [5] pre- 
sented a powerful new method , the Lyapunov weighted 
dynamics (LWD), which is not limited to low dimensions 
or restricted to a small family of problems. In this evolu- 
tionary approach, a swarm of walkers progress according 
to the rules of the underlying dynamics. The walkers 
proliferate or die depending on the degree of chaos en- 
countered by the system along a particular trajectory 
and, after many generations, only walkers on trajectories 
with the desired stability properties survive. Tailleur and 
Kurchan have demonstrated that their method is capable 
of finding even very small stability regions in systems of 
many degrees of freedom. 

Inspired by the work of Tailleur and Kurchan, we in- 
troduce here a general and efficient algorithm for find- 
ing trajectories with atypical stability properties, which 
is equally applicably to stochastic and deterministic dy- 
namics. The central notion of our approach lies in the 



definition of a Lyapunov weighted path ensemble, in 
which the statistical weight of trajectories explicitly de- 
pends on a measure for the chaoticity of the underlying 
dynamics. The degree with which particularly chaotic 
trajectories are favored or disfavored depends on the 
value of a parameter that can be viewed as conjugate 
to the measure of chaos, for instance the Lyapunov ex- 
ponent. This ensemble of trajectories is then sampled 
using techniques borrowed from transition path sampling 
(TPS), a method originally developed to study rare tran- 
sitions between long-lived stable states in complex molec- 
ular systems Note that a related combination of 
transition path sampling with a Lyapun ov weighted ac- 
tion has been suggested before [Ifl, By construc- 
tion, the transition path sampling procedure generates 
a Markov chain of trajectories distributed according to 
the chosen bias function. Following the terminology of 
Tailleur and Kurchan Q , we call this approach Lyapunov 
weighted path sampling. Similar techniques were recently 
used by Chandler and collaborators to sample an ensem- 
ble of trajectories weighted by an order parameter de- 
scribing the mobility of particles in a system undergoing 
the glass transition 
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The remainder of this article is organized as follows. 
In Sec. |TT]we briefly introduce the concept of Lyapunov 
instability and discuss ways to probe chaotic dynamics 
along individual trajectories. In Sec. IIIII we deflne the 
Lyapunov weighted path ensemble and in Sec. IIVI we 
explain how it can be sampled with transition path sam- 
pling algorithms. In Sec. |V]this approach is then applied 
to detect particularly stable and unstable trajectories in 
various systems. Some conclusions are provided in Sec. 
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II. LYAPUNOV INSTABLITY 
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Chaotic dynamical systems are characterized by a 
strong sensitivity to small changes in the initial condi- 
tions. To quantify this concept, consider a dynamical 
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system, described by a set of N coupled first order ordi- 
nary differential equations, 



(1) 



where x denotes a point in the iV-dimensional phase 
space. The time evolution of a infinitesimally small de- 
viation Sx separating two close-by trajectories is then 
governed by the linearized equations of motion 



Sx{t) = D{x)Sx{t), 



(2) 



where D{x) — dF/dx is the Jacobi matrix of the sys- 
tem evaluated at x. In a chaotic system, two points in 
phase space, initially separated by 6x{0) at time t = 0, 
will lead to trajectories that, on the average, separate 
exponentially in time, ~ |(5a;(0)| exp(At). Here, 

the vertical lines denote the Euclidean norm of a vector. 
The coefficient A, the long-time averaged growth rate of 
infinitesimally small displacements defined as 
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is called the Lyapunov exponent of the system. A positive 
Lyapunov exponent corresponds to exponential growth of 
an initially small perturbation and implies information 
loss and strong sensitivity to initial conditions, the defin- 
ing feature of chaotic dynamics. (It is possible to define 
whole spectra of Lyapunov exponents characterizing the 
growth rates of small perturbation in different directions 
of phase space [13|, In this article, however, we will 
consider only the largest Lyapunov exponent defined in 
Equ. ©.) 

Since Lyapunov exponents are defined as long-time av- 
erages, in an ergodic system every initial condition will 
yield the same A. For finite periods of time however, 
trajectories can display very different chaotic properties. 
An example are trajectories in "sticky regions" in the 
phase space of the standard map, in which trajectories 
can spend a long time before escaping away into more 
chaotic regions [isl [l6j . Other parts of phase space may 
be filled with regular periodic orbits that have a van- 
ishing Lyapunov exponents and are dynamically discon- 
nected with the chaotic surroundings. To describe such 
behavior of individual trajectories of finite length it is 
convenient to consider so-called finite time Lyapunov ex- 
ponents Xf{xo^t) that depend on the initial condition xq 
and on the temporal trajectory length 



(4) 



Such finite time Lyapunov exponents can be used to 
quantify the chaoticity of finite length trajectories. 

One difficulty occurring in the definition of the finite 
time Lyapunov exponent of Equ. (|4]), however, is that 
Xf also depends on the initial orientation of the dis- 
placement vector Sx(0). This vector reorients into the 
direction of fastest growth, but, depending on the degree 



of chaos prevalent in the respective phase space region, 
this reorientation may take a time long with respect to 
the trajectory length t. In fact, the time it takes to 
turn the displacement vector into the direction of the 
fastest growth is inversely proportional to the difference 
of the two largest Lyapunov exponents cx l/(Ai — A2) 
Here, Ai and A2 are the largest and the second 
largest Lyapunov exponent, respectively. This ambigu- 
ity in the definition of the finite time Lyapunov exponent 
can be avoided by integrating the equations of motion 
backwards for a time t_ longer than starting from the 
inititial condition xq . If the equations of motion for the 
displacement vector Sx are then integrated forward start- 
ing from X-t- with an arbitrary orientation of the dis- 
placement vector and following the reference trajectory, 
the displacement vector has oriented into the direction of 
fastest growth when xq is reached. Then, the displace- 
ment vector has a well defined orientation at i = and 
the definition of the finite time Lyapunov exponent for 
the trajectory from xq to Xt is unique. Since the reorien- 
tation time Tr may be large, this procedure can require 
the computation of long additional trajectory segments 
causing large computational costs. 

An alternative way to quantify the chaoticity of in- 
dividual finite length trajectories consists in determining 
the Relative Lyapunov Indicator (RLI) (l8l.lT9j. originally 
introduced to detect chaotic dynamics in planetary sys- 
tems. This measure has been proven particularly useful 
to distinguish regular from chaotic trajectories in sys- 
tems that are only weakly chaotic. The main idea of the 
RLI is to exploit the fact that in the chaotic regions of 
phase space, finite time Lyapunov exponents vary discon- 
tinuously as a function of the initial condition, i.e., adja- 
cent points can have very different local expansion rates 
[20L [21I . The RLI is defined as the magnitude of the dif- 
ference AA(a;o, t) between the finite time Lyapunov expo- 
nents of two trajectories separated by a small but finite 
amount Axq, 



AA(a;o,t) = \\{xo + /\xo,t) - X{xo,t)\ 



(5) 



It has been shown that the RLI is insensitive both to the 
separation Aiq as well as to the initial infinitesimal phase 
space displacement (5a; (0) used to calculate the finite time 
Lyapunov exponents A(a;o, t). To reduce fluctuations, one 
can average the RLI over time, 



t/At 



R{xo,t) = AX{xQ,At) 



^AA(xo,zAi), (6) 



where At is the time step used for the numerical integra- 
tion of the equations of motion. In the following we will 
use this smoothed version of the RLI to characterize the 
chaoticity of individual trajectories in various dynamical 
systems. 
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III. LYAPUNOV WEIGHTED TRAJECTORY 
ENSEMBLE 

As outlined in the Introduction, the goal of the work 
presented in this paper is to identify trajectories that 
have particular stability properties. For instance, one 
might be interested in locating regions of phase space 
that are populated by regular trajectories or in finding 
those pathways that are most chaotic. To do that, we 
start by defining an ensemble of trajectories including an 
additional weight that favors trajectories with the desired 
chaoticity properties. We assume that the dynamics gen- 
erates a stationary distribution p{x). Since we consider 
only deterministic systems here, we represent trajecto- 
ries of length t by their initial phase space point xq and 
define the Lyapunov weighted path ensemble (LWPE) as 

PL(xo,i) = ^PK)e"*«(^-*', (7) 

where the factor Q — f (ia;op(a;o)e"'^^^"'*'' normalizes the 
distribution. Here, the parameter a, which can be viewed 
as conjugate to the chaoticity indicator R{xo, t), controls 
how strongly the weight of initial condition xq is changed 
according to the chaoticity of the trajectory evolving out 
of Xq. An analogous ensemble of trajectories can be easily 
defined using finite time Lyapunov exponents if they can 
be calculated accurately (simply replacing R{xo,t) with 
Xf{xo, t) in the above equation). Large positive values of 
a favor strongly chaotic trajectories with a large chaotic- 
ity indicator. For negative as, on the other hand, weakly 
chaotic or regular trajectories with a small chaoticity in- 
dicator are given a larger weight in the ensemble. Note 
that a similar ensemble of trajectories can easily be con- 
structed also for stochastic dynamics and appropriately 
defined Lyapunov exponents. 



IV. SAMPLING TRAJECTORIES WITH TPS 

We sample the Lyapunov weighted path ensemble 
with a technique borrowed from transition path sam- 
pling d, H, Q, a method originally developed to simu- 
late rare but important transitions between long-lived 
stable states as they occur, for instance, in protein fold- 
ing, chemical reactions and first order phase transitions. 
In a transition path sampling simulation a biased ran- 
dom walk is carried out in the space of trajectories in a 
way such that trajectories are sampled according to their 
weight in the desired ensemble. This is accomplished us- 
ing a Monte Carlo procedure in which a trial trajectory is 
generated from the current trajectory and then accepted 
according to the Metropolis rule. Iterating this basic 
step, a set of trajectories with the correct probability is 
generated. A particularly efficient way to generate trial 
pathways is the so called shooting algorithm 22, 23], also 
used in the present work. In this approach, a new tra- 
jectory is generated from an old one by first randomly 



selecting a point on the old trajectory, and then inte- 
grating the equations of motion starting with perturbed 
momenta. The magnitude of the perturbation of the mo- 
menta controls how different the new trajectory is from 
the old one and, therefore, also controls the average ac- 
ceptance probability of the Monte Carlo procedure. 

In more detail, the path sampling procedure of the 
Lyapunov weighted path ensemble is carried out in the 
following way. The first trajectory is created by integrat- 
ing the equations of motion starting from an arbitrary 
initial condition xq. From this trajectory one then selects 
a point at random. At this so-called shooting point, the 
momenta are slightly changed by addition of a small per- 
turbation drawn from a Gaussian distribution. The new 
trajectory is then obtained by integrating the equations 
of motion forward to time t and backward to time 0. The 
new trajectory with initial condition x^'q^ and chaoticity 

indicator R{xQ^\t) is accepted with probability 

I p(4 ) J 

where Xg"'' and R{xl^\ t) denote the initial point and the 
chaoticity indicator of the old trajectory, respectively. If 
the new trajectory is accepted, the procedure will be re- 
peated with this trajectory. Otherwise the old trajectory 
is kept as the current one. The acceptance probability of 
Equ. (Ill) is derived from the detailed balance condition 
and guarantees that trajectories are harvested according 
to the Lyapunov weighted path ensemble. Note that to 
obtain Equ. ([8] ) we have assumed a momentum pertur- 
bation that leads to a symmetric generation probability. 
If this is not the case, an appropriate factor must be 
taken into account in the acceptance probability. The 
shifting algorithm of transition path sampling [2^, [2^ 
can be adapted in a similar way to sample the Lyapunov 
weighted path ensemble. In the next section we discuss 
the application of this method to various chaotic dynam- 
ical systems. 



V. RESULTS 

In this Section we use the method outlined above to 
identify particularly stable and unstable trajectories in 
various dynamical systems with dimensionality ranging 
from two to several hundred. For better comparison, we 
mostly follow Rcf. [5] in our choice of examples. 



A. Standard map 

The standard map is a representation of the dynamics 
of a free rotor kicked at regular intervals with an impul- 
sive force in a given direction [l|, [23| . For a kicking period 
of 1 and a kick strength of fc, the standard map is given 
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by 



Wn+i =0Jn- — sin(27r(p„), 

ZTT 



(9) 



where a;„ and are the angle and the angular momen- 
tum of the rotor immediately after the n-th kick, respec- 
tively. Due to the periodicity of the motion, both variable 
Lo and Lp are considered on a torus (by taking lo and Lp 
modulo 1). The standard map is area preserving and, 
depending on the value of k, displays different degrees of 
chaos. For vanishing k, the dynamics reduces to the mo- 
tion of a free rotor and the system is integrable. Accord- 
ingly, the largest Lyapunov exponent vanishes. As the 
parameter k is turned on, some chaos develops in par- 
ticular regions of the two-dimensional phase space, while 
other regions originating from KAM-tori remain regular 
resulting in a phase space structure consisting of islands 
of stability embedded in a chaotic sea P, H. As k IS 
increased, the stable regions shrink and for large values 
of k only a very small fraction of initial conditions lead 
to regular trajectories. In the following we will sample 
the Lyapunov weighted path ensemble to find these rare 
regular trajectories. 

To sample trajectories of the standard map we use 
the standard shooting algorithm, in which a new tra- 
jectory is obtained from the current one by first selecting 
a phase space point on the current trajectory. Then, the 
selected point is slightly perturbed by addition of a ran- 
dom displacement drawn from a Gaussian distribution 
with width ctq to both Lp and u. Starting form this per- 
turbed initial condition, the new trajectory is obtained 
by carrying out an appropriate number of iterations of 
the standard map in forward direction and of the inverse 
map in backward direction. Finally, the new trajectory 
is accepted with the acceptance probability of Equ. (HJ. 



Using the shooting algorithm we have generated 10^ 
trajectories of length n = 10* with o-q = 0.05 for a 
kicking strength of fc = 7.7. In this case, the chaotic 
sea covers almost the entire phase space. The magni- 
tude of initial point deviation for the RLI calculation 
was |Aa;o| = 10"-^^ and the probability density /o(a;o) of 
initial conditions was assumed to be uniform. To find the 
rare regular trajectories a negative value of the control 
parameter a = —4 was used. The right choice of this pa- 
rameters is crucial as it controls how phase space is sam- 
pled in a way that is analogous to the effect of the inverse 
temperature in a regular Monte Carlo simulation. A very 
large value of a favors regular trajectories very strongly 
such that the simulation gets trapped very easily in local 
minima of the RLI and the most regular trajectories are 
not found. A small value of a, on the other hand, does 
not generate a sufficiently strong bias towards the atypi- 
cal regular trajectories such they may not be found in this 
case either. Using these parameters, several rare regular 
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FIG. 1. (a) Phase space plot of the standard map for k = 
7.7. The arrows indicate the locations of the stability islands, 
which are not visible at the scale of the figure, (b) Enlarged 
phase space portraits of the regular trajectories found by the 
Lyapunov weighted path sampling algorithm. 



trajectories were found in the simulation as displayed in 
Fig. [H These regular trajectories correspond to the four 
stable islands found also in Ref. Q • For the stable orbits 
the RLI values ranged from 10~^* to 10~^^ compared to 
values of 10"'' for typical trajectories in the chaotic sea. 
As can be inferred from Fig. [TJ the regular islands cover 
only a very small part of phase space. The fact that reg- 
ular trajectories are nevertheless found in the simulation 
indicates that the landscape of the chaoticity indicator 
must have a global funnel-like shape that attracts the 
simulation towards the regions of highest regularity. 



B. Spring pendulum 

We next consider the spring pendulum evolving ac- 
cording to Hamilton's equation of motion. This two- 
dimensional system consists of a point of mass m exposed 
to a constant force of magnitude g in negative y-direction 
attached to a fixed pivot by a harmonic spring such that 
both angle and length of the pendulum can change in 
time. The Hamiltonian of the spring pendulum is given 
by 



(10) 
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FIG. 2. Very stable and unstable trajectories of the spring 
pendulum traced out in configuration space. In both cases 
the energy is _E = 2. 

where r ~ \J ^ is the distance of the mass point 
to the pivot, = + is the squared magnitude of 
the momentum, k is the spring constant and K is the 
equilibrium length of the spring [2^ [25|. In all follow- 
ing calculations, the mass m, the equilibrium length i? 
and the force constant k are set to unity and the force 
strength to g = 2. 

For this model, we have carried out a Lyapunov 
weighted path ensemble simulation for two different val- 
ues of the parameter a. In both cases, the total en- 
ergy was E = 2 and the equations of motion are inte- 
grated with a time step of Ai = 10~^ using the sym- 
plectic time-reversible Forest-Ruth algorithm of fourth 
order [2^ yielding good energy conservation. The shoot- 
ing algorithm was carried out in the usual way by adding 
perturbations drawn from a Gaussian distribution to the 
momenta at the shooting point. Results of these simula- 
tions are displayed in Fig. [2j The very stable trajectory 
of length T — 100 was obtained after 4000 iterations of 
the path sampling algorithm with a = — 5x10^^ favoring 
regular trajectories. In this calculation the magnitude of 
the shooting displacement was ug = 0.5. This trajec- 
tory has an RLI of ^ 2 x 10~^^, a value which is about 
four orders of magnitude smaller than that of a typical 
trajectory. 

The unstable trajectory of Fig. [2]has been found after 
3000 iterations carried out with a — 10^ and a shooting 
displacement of magnitude cjq — 5 x 10""*. This trajec- 
tory has an RLI of~ 6x lO^"', which is about 6 orders of 
magnitude larger than that of typical trajectories. The 
qualitative difference between the two trajectories of Fig. 
[2] is remarkable. While the stable trajectory has much 
amplitude in the angular degree of freedom, the unsta- 
ble one displays pronounced stretching movements with 
smaller angular oscillations. 

C. Fermi- Pasta-Ulam chain 

As an example of a system with higher dimensional- 
ity, we search for particularly stable trajectories of the 
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FIG. 3. Particle positions of a FPU-chain with A'^ = 32 as 
a function of time along part of a trajectory sampled by the 
Lyapunov weighted path sampling algorithm. Particle posi- 
tions have been displaced vertically for better visibility. The 
trajectory was obtained after a couple of thousand iterations 
with Q = — 5 x 10* favoring particularly stable trajectories. 



Fermi-Pasta-Ulam (FPU) chain, originally concocted to 
study the thermalization of oscillatory modes in solids 
(27I l28j . This model consists of a one-dimensional chain 
of N point particles with mass m located at positions Xi . 
The particles are coupled by harmonic springs to which 
a weak anharmonic part in form of a quartic potential is 
added leading to the Hamiltonian 

i i 

Here, k is the spring constant of the harmonic spring and 
/3 is a parameter controlling the strength of the quartic 
potential. Below, we use units in which m = 1, fc = 1 and 
consider the case /? = 0.1 with fixed boundary conditions. 

Depending on initial conditions, the FPU-chain dis- 
plays different types of motion. Starting from a state in 
which all the energy is concentrated in a high frequency 
mode, the system evolves into a so-called breather mov- 
ing chaotically with the energy strongly localized in space 
before equilibration eventually sets in on very long time 
scales [29|. Other initial conditions lead to only weakly 
chaotic solitonic modes, in which a kink moves through 
the system with constant speed and a preserved shape 
of its sharp front [28;, Here, we use our algorithm 
to find these weakly chaotic solitons by using a param- 
eter a — —5 X 10^ highly favoring regular trajectories. 
Note that this strongly negative value of a was neces- 
sary since chaos in the FPU-chain in this regime is very 
weak and it is difficult to distinguish between more and 
less chaotic trajectories. Particle traces for a system of 

= 32 particles at a total energy oi E = 32 are shown in 
Fig. [3) The equations of motion were integrated with a 
time step of At — 0.05 for trajectories of length r = 5000. 
The trajectory shown in Fig. [3] was obtained after a cou- 
ple of thousand iterations of the path sampling scheme 
with shooting displacements in momentum space of mag- 
nitude CTG = 0-5. Kinks wandering through the chain at 
constant velocity and bouncing back and forth between 
the chain ends are clearly visible indicating a solitonic 
mode of motion. 
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D. Double well system 

To test whether Lyapunov biased path samphng can 
be used to find reactive trajectories in systems with mul- 
tiple stable states separated by barriers, we have studied 
a simple double well system in one dimension. The dy- 
namics of the system is governed by the Hamiltonian 



Zm 



(12) 



corresponding to a particle of mass m moving on a po- 
tential energy surface with minima at x ~ ±l,y — 
and a maximum at the origin. The potential energy 
barrier separating the two stable states has a height of 
k. We imagine that the system is in contact with a 
heat bath with temperature T, such that the distribu- 
tion of initial conditions of the trajectories is canonical, 
p{x,p) oc exp{— Here, (3 — l/k^T is the re- 
ciprocal temperature and fee is the Boltzmann constant. 
Thus, all initial conditions of the system are in princi- 
ple accessible in this ensemble, albeit with different sta- 
tistical weight. This system is integrable and hence its 
Lyapunov exponent vanishes. Nevertheless, finite length 
trajectories diverge strongly near the barrier top. 

Although a canonical distribution of initial conditions 
implies a coupling of the dynamics to the degrees of free- 
dom of the heat bath, we assume that this coupling is so 
weak that it does not affect the dynamics of the system on 
the time scale of the length r of the trajectories. Hence, 
each individual trajectory evolves at constant total en- 
ergy. Under these conditions, the acceptance probability 
of the path sampling procedure is given by 
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(13) 



where H{-b) and H{o) are the energies of the new and the 
old trajectories, respectively, and i?(n, r) and i?(o, r) are 
the respective RLIs. The above acceptance probability 
takes into account possible energy changes due to the 
momentum perturbation applied at the shooting point. 

Since we are interested in reactive trajectories that 
cross the energetic barrier between the stable states and 
the phase space structures near saddle points which lie 
at origin of chaos jl| , we sample the Lyapunov weighted 
path ensemble procedure with a large value of a = 5 x 10^ 
strongly favoring chaotic trajectories. The equations of 
motion are integrated with a time step of At — 10^^ 
and trajectories have a total temporal length of r = 10^. 
Units where chosen such that m = 1 and k = \ and the 
temperature was set to /3 = 1. The magnitude of the 
shooting displacement to the momenta was <jq = 0.05. 
Starting from a trajectory of energy E — 0.05 oscillating 
about the bottom of one well, after a few hundred path 
sampling steps reactive trajectories connecting the two 
wells were obtained. The phase space plot of a reactive 
trajectory with total energy E = 1.0000019 (just slightly 
above the barrier height oi E = 1.0) crossing the barrier 
in close proximity of the saddle point is shown in Fig. 2] 




FIG. 4. Phase space plot of a particularly chaotic trajectory 
crossing the saddle point in the double well system. Note that 
at the origin there is a gap between the upper and the lower 
branch of the periodic trajectory, but the gap is too small to 
be visible at the scale of the figure. 



E. Double well dimer in a solvent 

As exemplified by the results of the previous section, 
the dynamical character of pathways crossing barriers 
in the vicinity of saddle points in the potential energy 
surface strongly differ from that of trajectories fluctuat- 
ing about minima Such pathways connecting stable 
states are, for instance, relevant in the context of ac- 
tivated chemical reactions and there is intense interest 
in computational methods for finding such rare barrier 
crossing pathways along which chemical reactions occurs 
Q . In this section we will study if the Lyapunov weighted 
path sampling method can be used to identify reactive 
trajectories based on their chaoticity. In particular, we 
will address the question wether such an approach can be 
successful for reactions occurring in solution, where the 
chaoticity of the reactive subsystem may be overshad- 
owed by that of the solvent. 

We study this issue using a simple two-state model 
dimer embedded in a soft-sphere solvent [3l|, [s^ . This 
three-dimensional model consists of N particles of mass 
m evolving according to Hamilton's equations of motion 
in a cubic box of volume V and periodic boundary con- 
ditions. All particle interact pairwise via the purely re- 
pulsive Weeks-Chandler- Andersen (WCA) potential [ss} . 



VwcA{r) 




e for r < Tc 
for r > Tr 



(14) 



Here, r is the interparticle distance, cr is the interaction 
radius, e is the strength of the potential, and Tc = 
is the cutoff radius. In addition, the two particles forming 
the dimer are bonded through the double well potential 



Kiw(r) 
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(15) 



where h and w are constants determining the height and 
the width of the barrier, respectively. The distance be- 
tween the two minima of the double well potential is 
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FIG. 5. Interatomic distance r of the double well dimer as a 
function of time t along reactive trajectories found by the Lya- 
punov weighted path sampling procedure for particle numbers 
iV = 22 (top), iV = 32 (center), and iV = 108 (bottom). In 
all three trajectories a transition from the contact to the ex- 
tended state occurs as evidenced by the sudden increase of r 
towards the end of the trajectory. 



Iw. For low temperatures, i.e., for h ^ fceT, the dimer 
mainly exists in two states. In one state, the contact 
state, the interatomic distance of the dimer fluctuates 
about Tc- In the other state, the extended state, the 
dimer has a bond length of about Tc + 2i/;. Thermally 
activated transitions between these two states occur very 
rarely and are separated by long permanence times in the 
wells. 

In all our simulations we use reduced units in which 
CT = 1, e = 1, and m — 1. We study this system for 
particle numbers iV = 22, = 32 and N = 108 at a 
density of p = N/V — 0.5 and at a total energy per 
particle of E/N = 1.0 corresponding to a temperature 
of about T = 0.55. The dimer barrier height was set 
to h — 6.0 and its width to w — 0.3. Trajectories of 
total length r = 10 were integrated with a time step of 
At — 10^'^. In the path sampling procedure, different pa- 
rameters a and different magnitudes ctq of the shooting 
displacements were used. While for the smallest system 
a = 170 and ctg = 0.2 were employed, we chose a ~ 160 
and CTG = 0.25 for the system of intermediate size and 
a = 600 and ctg = 0.05 for the largest system. Starting 
from non-reactive trajectories, the Lyapunov weighted 
path sampling algorithm carried out with these parame- 
ters succeeded in finding reactive trajectory. Examples of 
the interatomic distance as a function of time are shown 
in Fig. [S] along reactive trajectories for various system 
sizes. While for all system sizes studied here the Lya- 
punov weighted path sampling simulation converged to- 
wards reactive trajectories, the number of iterations be- 
fore observing the first reactive path increased with sys- 
tem size. While for particle number N = 22 the first reac- 
tive event occurred after about 400 path sampling steps. 



the first reactive trajectory was found after 1000 itera- 
tions for N = i2 and after 2900 iterations for N = 108. 
These results indicate, that the pronounced chaoticity 
of barrier crossing trajectories can be used to identify 
reactive pathways even for systems with hundreds of de- 
grees of freedom. Note, however, that not all reactive 
trajectories are strongly chaotic. During our simulations 
it happened repeatedly that reactive pathways were re- 
jected, because their relative Lyapunov indicators were 
too small. It could be that along these trajectories the 
chaoticity associated with the saddle point crossing was 
compensated by a particularly stable dynamics of the 
solvent. 



VI. CONCLUSIONS 

In this paper, we have presented a fiexible numerical 
method to find particularly chaotic or regular trajecto- 
ries in dynamical systems. The basic idea of the method, 
inspired by the work of Tailleur and Kurchan Q and 
called Lyapunov weighted path sampling, is to first de- 
fine an ensemble of trajectories weighted by a measure of 
their chaoticity, for instance their finite time Lyapunov 
exponent. In this trajectory ensemble a parameter, which 
can be viewed as conjugate to the Lyapunov exponent, 
can be tuned to favor either very chaotic or very regular 
trajectories. The trajectory ensemble is then sampled 
with methods adopted from transition path sampling. 
Other chaoticity indicators besides finite time Lyapunov 
exponents can be easily integrated into the algorithm as 
well. Since the calculation of finite time Lyapunov expo- 
nents can be computationally demanding, we have, for 
instance, used relative Lyapunov indicators (RLI) to bias 
trajectories according to their level of chaos. These in- 
dicators are particularly sensitive and are capable of dis- 
tinguishing weakly chaotic trajectories from regular ones. 
While in this paper we have used Lyapunov weighted 
path sampling to study only systems evolving determin- 
istically, the method can be applied as easily to stochastic 
dynamics provided an appropriate chaoticity indicator is 
available. 

The complexity of the examples studied here ranges 
from a simple one-dimensional double well system to the 
FPU-model and a bistable dimer in a solvent with hun- 
dreds of degrees of freedom. In all cases, the Lyapunov 
weighted path sampling algorithm successfully identified 
trajectories with atypical chaoticity properties. While for 
the FPU-model we used Lyapunov weighted path sam- 
pling to find weakly chaotic solitonic modes of motion, we 
concentrated on highly chaotic trajectories for the dimer 
in solution. The results obtained for this simple model 
of a chemical reaction indicate that it is possible to use 
Lyapunov weighted path sampling to find rare reactive 
trajectories that pass through saddle points in the po- 
tential energy surface as they connect long-lived stable 
states with each other. Further studies will be neces- 
sary to clarify to which degree identifying such trajecto- 
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ries is made difficult by the chaos arising from degrees 
of freedom not directly coupled to the reaction (for in- 
stance solvent degrees of freedom) and possibly eclipsing 
the dynamical instability of the reactive subsystem. It 
will also be interesting to investigate whether chaotic- 
ity indicators such as the maximum Lyapunov exponent, 
the Kolmogorov-Sinai entropy or the relative Lyapunov 
indicators used in the present study correlate with the 
measures of mobility used by Chandler and collabora- 
tors to link the glass transition with a first order phase 
transition in trajectory space [l2|- In their work, these 
authors started from the equilibrium distribution of path- 
ways and added to it a bias that favors trajectories with 
low dynamical activity. Chandler an coworkers demon- 
strated numerically that this transition displays all the 
features of a first-order transition occurring in trajectory 
space. It would be interesting to study if an analogous 



bias based on chaoticity indicators also leads to an equiv- 
alent first order transition in path space. In such research 
it may be fruitful to combine Lyapunov weighted path 
sampling with advanced equilibrium simulation methods 
such as umbrella sampli ng j34i| . metadynamics [35j], or 
parallel replica sampling [36| directly acting on chaotic- 
ity indicators. 
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